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. We present an efficient scheme for representing many-body wavefunctions in quantum Monte 

, ' Carlo (QMC) calculations. The scheme is based on B-splines (blip functions), which consist of 

^ , localized cubic splines centred on the points of a regular grid. We show that blip functions are 

• unbiased, systematically improvable, and conveniently obtained from any standard plane-waves 

' density functional theory (PW-DFT) code, and therefore provide a convenient and natural interface 

between PW-DFT and QMC calculations. We present tests on a 16-atom system of Si in the /3-tin 
' ^ I ' ■ structure, and on 2- and 8- atoms systems of MgO in the NaCl structure. We show that already 

O ' with such small systems the speed-up of blip functions with respect to plane-waves is between one 

, and two order of magnitudes, without compromising the accuracy. 

r^H ' The quantum Monte Carlo (QMC) technique [1] is becoming ever more important in the study of condensed matter, 
, with recent applications including the reconstruction of semiconductor surfaces [2], the energetics of point defects in 
^ insulators [3], optical excitations in nanostructures [4], and the energetics of organic molecules [5]. Although its 
demands on computer power are much greater than those of widely used techniques such as density functional theory 
(DFT), its accuracy is also much greater for most systems. With QMC now being applied to large complex systems 
containing hundreds of atoms, a major issue is the basis set used to represent the many-body wavefunctions. We 
propose and test here a basis set that has many of the properties of plane waves, currently widely used in density 
^ , functional theory (DFT) calculations, but with the advantage of being localized. 
Q ■ In QMC, the trial many-body wavefunction ^'T(ri, . . . v^) consists of a Slater determinant D — or more generally a 
^ linear combination of Slater determinants — of single-electron orbitals ■i/'„(ri) multiplied by a parameterized Jastrow 
correlation factor J(ri, . . .r^)- In variational Monte Carlo (VMC), J is 'optimized' by varying its parameters so as 

■ to reduce the variance of the 'local energy' 5*^^ (^H'^t^, where H is the many-electron Hamiltonian. Since VMC by 

itself is not usually accurate enough, the optimized produced by VMC is used in diffusion Monte Carlo (DMC), 
CO , which achieves the exact ground state within the fixed nodal structure imposed by the Slater determinant D. At 

■ each QMC step it is necessary to evaluate \I't(i"i, . . .tat) in each of the replicas (QMC "walkers"), which involves 
the evaluation of the single electron orbitals ^„(ri). A crucial issue in the efficiency of the calculations is therefore 
the representation of ipn(j^i)- One common approach is using plane waves (PW). The big advantages of PW are that 

I their accuracy is systematically improvable (by increasing the PW cutoff) and they are unbiased. Moreover, many 
DFT codes are written in terms of PW, so the technology is highly developed and easily accessible. However, PW 
d ' are not well suited for QMC calculations because for the evaluation of each ipn{j^i) one has to sum over all PW in the 
system. Since this has to be done for M orbitals and N electrons, with M proportional to N, the cost of evaluating 
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I ' the many-body wavefunction involves 0{N^) operations, with a prefactor depending on PW cutoff, which can be very 



large for 'hard' systems, like MgO. The storage required for a PW representation is proportional to iV^. 
Q ' This problem with PW can be overcome by using localized basis sets. One possibility is to use gaussians, but the 
O ■ drawback is that they are biased and generally difficult to improve systematically. An option which combines the 
best of both worlds, is to use a B-spline basis (blip functions), already proposed for 0{N) DFT calculations [6]. Here 
we propose and test the use of blip functions in QMC calculations. We will show that blip functions share all the 
/\ ' advantages of PW, i.e. are systematically improvable and unbiased. They are also closely connected with PW, and 
^ I can therefore be obtained from PW-DFT codes. However, they are localized, therefore the evaluation of each orbital 
ipni^i) has a cost which is independent of the size of the system and indeed of blip-grid spacing (connected to the PW 
cutoff). The storage required for blip functions is not much worse than PW and has the same 0{N^) scaling. 

As described in detail elsewhere [6] , the blip functions consist of localized cubic splines centred on the points of a 
regular grid, each function being non-zero only inside a region extending two grid spacings in each direction from its 
centre. For a cubic grid spacing a, the blip fmiction Os(r) centred on the grid point at position = {Xg, Yg, Zs) is 
given by: 

e,(r) = ^{{x - Xg)/a)cj){{y - i;)/a)0((z - Z,)/a), (1) 

where (^(^) is: 
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= J(2-i?if i<iei<2 

= |e|>2. (2) 

The function and its first two derivatives are continuous, discontinuities appear only in the third derivative, and all 
higher derivatives are zero. Each single particle orbital is then represented as: 

V'nW = X^ftnsQsW- (3) 

For any position r, there are only 64 non-zero blip functions, whatever the nature and size of the system, so that the 
number of operations to compute V'n(r) is the same for any material. 

The close relationship between B-splines and PW has been discussed elsewhere [6]. In the PW representation, the 
single particle orbitals are given by: 

V'„(r)=^c„ke*-- (4) 

k 

where the wavevectors k go over the reciprocal lattice vectors of the superlattice, with k less than the PW cutoff kmax- 
The relationship between the PW coefficients c„k and the blip coefficients a„s can be understood by considering blip 
waves Xk(r) defined by: 

Xk(r) = ^e^''«-e,(r). (5) 

S 

For small k, the Xk(r) are essentially identical to plane waves exp(k • r), apart from a k-dependent factor 7k: 

e^'''-^7kXk(r). (6) 

The factor 7k is the Fourier transform of a single blip function 0(r) and is given by 7k = 'Jkxlkylk^, where k = 
{kx,ky,kz) and 

3 

= — (3-4cosA;-|-cos2fc). (7) 

At larger k, the Xk(i") differ significantly from exp(k • r), as they must, because XkCr) is periodic in k space: the 
number of independent Xk(r) functions is equal to the number of sites on the blip grid. 

There is a "natural" choice of blip grid spacing a, given by a = Tr/Zcmax- With this choice, the region k ~ kmax 
where blip-waves and plane-waves differ most is the region where the plane waves coefficients c„k are very small. 
However, the precision with which blip-waves reproduce plane-waves in the region k < kmax can always be improved 
by refining the blip grid. 

The procedure to obtain the blip coefficients a„s from the plane- wave coefficients of orbitals Vn(r) obtained from 
a DFT calculation is straightforward. For the relationship between blip- waves and plane- waves (see eqn. 3, 4 and 7), 
it follows that 



,k7ke'''■"^ (8) 

The coefficients 

CLfis can therefore be evaluated using fast Fourier transform routines. 

We have implemented blip functions in the appropriately modified CASINO code [8] . To test the implementation we 
present now three cases in which we compare the energy and the standard deviation in VMC and DMC calculations 
performed using PW or blip-functions representations of the single- particle orbitals. Calculations with blip functions 
are presented for two values of the grid spacing, the natural grid spacing a = Tr/fcmax and a two times finer grid 
obtained with a = 7r/2fci„ax- Results are reported in Table I. All calculations have been performed at the F point. 

The first case is a 16-atom cell of silicon in the /3-tin structure. The single-particle orbitals have been obtained 
using the pwscf code [9], with Hartree-Fock pseudo-potentials (p-channel chosen to be the local part) and PW cutoff 
energy of 15 Ry. VMC calculations are reported for 3.2 x 10^ steps of length 1 a.u. in all three cases. No Jastrow 
factor has been used for these VMC calculations. DMC calculations have been performed using 320 walkers for 10100, 
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12700 and 10100 steps of length 0.03 a.u. for PW and blip functions calculations with the coarse and the fine grid 
spacing respectively. Diffusion to the ground state is already achieved after ^ 100 steps. In VMC the natural grid is 
not dense enough for this system, with the largest difference being in the kinetic energy (of Ri 0.06 eV/atom). The 
standard deviation on the energy is also slightly larger. However, with the fine grid the blip functions results agree 
identically with the PW ones within a statistical error of only a few meV/atom. In Table I we also report the time 
taken to perform one VMC step on an Origin 3000 machine. Already for this small system, with such a modest PW 
cutoff, the speed-up with blip functions is almost a factor of 6. The timings between the two blip- function calculations 
should in principle be identical, the small difference between the two is probably due to the larger sparsity in memory 
of the blip coefficients a„i for the case with a finer grid, and we found that this is machine dependent. For DMC the 
computational speed-up is more that a factor of 10, and the energy is already correct with the natural grid, which 
means that the nodal surface is essentially the same as the PW one already with the natural grid. 

The second test we performed was a perfect crystal of MgO in its zero pressure NaCl structure. The unit cell in 
this case contained only 2 atoms and had face-centred-cubic (f.c.c.) geometry. Single-particle orbitals were obtained 
again using the pwscf code, with Hartree-Fock pseudo-potentials (d-channel chosen as the local part for both Mg 
and O) and a PW cutoff of 200 Ry. No Jastrow factor has been used in these calculations. VMC calculations have 
been done with 3.36 x 10^, 1.6 x 10^ and 1.6 x 10* steps of length 0.3 a.u. for PW and the two blip functions cases. 
DMC calculations have been performed using 1600 walkers for 10.79 x 10'', 23.98 x 10'', 11.53 x 10"* steps of length 
0.005 a.u. for the three cases [7]. Diffusion to the ground state is achieved after the first few hundred steps. Similarly 
to the previous case, blip functions VMC energies and standard deviation agree identically with those obtained using 
PW for the dense blip grid, and the speed-up obtained with blip functions is more than a factor of 10. DMC energies 
are also in this case correct already when the coarse grid is used, but the variance is significantly improved when the 
fine grid is used. 

Finally, the third test consists of the same MgO crystal simulated in a simple cubic (s.c.) cell containing 8 atoms. 
Single-particle orbitals were obtained in analogy to the previous case, i.e. same pseudo-potential and same PW cutoff 

of 200 Ry. No Jastrow factor has been used. VMC calculations have been done with 0.32 x 10®, 1.6 x 10® and 
1.6 X 10® steps of length 0.3 a.u. for PW and the two blip functions cases. The important thing to notice in this case 
is the speed-up obtained with blip functions, which is over two order of magnitudes. We have not attempted DMC 
calculations as they would be impractical for the PW case. 

We note that despite we have chosen to use PW cutoffs of 15 Ry and 200 Ry for Si and MgO respectively, we found 
that by using much larger cutoff energy (typically 32 Ry for Si and 500 Ry for MgO) the variance of the energy can be 
further significantly improved. Of course, increasing the PW cutoff leads to a direct increase in the PW computational 
time, but has hardly any effect in the calculations which employ blip functions. We have also found that by using 
much larger PW cutoff the blip functions natural grid is already accurate enough, as expected. 

We have presented here a robust and efficient scheme based on B-splines to represent the trial wavefuntions in 
QMC calculations. We have shown that this scheme shares all the advantages of plane-waves, but offers a much 
better scaling behaviour with respect to the number of atoms in the system and the hardness of the pseudo-potentials 
used in the calculations. This scheme has been implemented in the CASINO [8] code, and we have presented tests on 
three different cases. The largest system considered here (in terms of number of plane waves) was an MgO crystal 
in the NaCl structure simulated with a s.c. imit cell containing 8 atoms. We have shown that already for this 
relativ(^ly small system the speed-up obtained using blip functions is over a factor of 100. Since B-splincs can easily 
be obtained from PW, they also provide a natural and convenient interface between QMC and PW-DFT codes. 
Moreover, this technique can be used in conjuction with "linear-scaling" techniques for QMC calculations, as reported 
elsewhere [10,11]. We conclude by noting that we are now attempting to calculate the formation energy of a Schottky 
defect in MgO using a cell containing 54 atoms. This calculation would be impossible to perform if we had to use 
PW (results will be reported elsewhere [12]). 
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PW 


Blips(a = 7r/fcmax) 


Blips(a = 7r/2A;max) 


Si /3-tin, 16 atoms 








VMC 








Ekin 
Eloc 
Enl 

Etot 
T (s/step) 

DMC 


43.864(3) 
15.057(3) 
1.533(3) 
-101.335(3) 4.50 
1.83 


43.924(3) 
15.063(3) 
1.525(3) 
-101.277(3) 4.74 
0.32 


43.862(3) 
15.058(3) 
1.535(3) 
-101.341(3) 4.55 
0.34 


Etot 
T (s/step) 


-105.713(3) 2.29 
2.28 


-105.711(4) 2.95 
0.21 


-105.715(4) 2.38 
0.25 


MgO-NaCl, 2 atoms, f.c.c. cell 








VMC 








Ekin 
Eloc 
Enl 

Etot 

T (s/step) 


199.449(24) 
-239.899(27) 
-26.906(12) 
-224.527(4) 28.7 
101 X 10"^ 


199.465(15) 
-239.861(15) 
-26.889(8) 
-224.465(3) 35.8 
8.3 X 10"^ 


199.418(15) 
-239.855(15) 
-26.902(8) 
-224.523(2) 28.3 
8.9 X 10"^ 


DMC 








-Etot 
T (s/step) 


-228.429(10) 22.1 
89 X 10"^ 


-228.433(7) 28.9 
7.1 X 10"^ 


-228.427(9) 22.3 
7.5 X 10"^ 


MgO-NaCI, 8 atoms, s.c. cell 








VMC 








Ekin 
Eloc 
Enl 

Etot 
T (s/stop) 


178.349(49) 
-225.191(50) 
-17.955(25) 
-227.677(8) 14 
7.8 


178.360(22) 

-225.128(24) 
-17.974(11) 
-227.648(4) 15 
5.6 X 10"^ 


178.369(22) 

-225.177(23) 
-17.976(11) 
-227.669(4) 14.5 
7.1 X 10"^ 



TABLE L Comparisons of the various components of the total energy (in eV/atom) and timings between VMC and DMC for 
a 16 atoms Si system in the /3-tin structure, and an MgO crystal in the NaCl structure. Standard deviation of the total energy 

(eV/atom) is also reported beside the total energy. The MgO crystal has been simulated using a 2 atoms face-centred-cubic 
cell and an 8 atoms simple cubic cell. PW calculations have been performed with a cutoff energy of 15 Ry for Si and 200 Ry 
for MgO. Blips calculations have boon performed using two different grid spacing: a = Tr/fcmax and a = 7r/2fcniax, where fcmax 
is the modulus of the largest PW vector. 
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